Appendix C: R script applied (modified from Silva et al. 2016)

Silva DP, Vilela B, Buzatto BA, Moczek AP, Hortal J (2016)
Contextualized niche shifts upon independent invasions by the dung beetle Onthophagus taurus. Biol Invasions 18:3137–3148

# ---------------- Code used in "Climatic niche evolution among native and introduced areas (as an archaeophyte and neophyte) of Lilium lancifolium (Liliaceae)" ----------------

# ---------------- Installation of packages ----------------

install.packages("BIOMOD",repos="http://R-Forge.R-project.org")
install.packages("ade4")
install.packages("adehabitat")
install.packages("sp")
install.packages("gam")
install.packages("MASS")
install.packages("mvtnorm")
install.packages("gbm")
install.packages("dismo")
install.packages("knitr")
install.packages("spThin")
install.packages("rgeos")
install.packages("maptools")
install.packages("raster")
install.packages("ecospat")
install.packages("spThin")
install.packages("rgdal")
install.packages("riverplot")

# ---------------- Importation of libraries ----------------

library(BIOMOD)
library(ade4)
library(adehabitat)
library(sp)
library(gam)
library(MASS)
library(mvtnorm)
library(gbm)
library(dismo)
library(knitr)
library(spThin)
library(rgeos)
library(maptools)
library(raster)
library(ecospat)
library(spThin)
library(rgdal)
library(riverplot)

# ---------------- Parameters to modify by the user ----------------

# Set working directory
main.directory <- "E:/Lilium lancifolium/TEXT/Supplementary/RCode"

#Defining groups (regions) to be tested and the colors for each region
n.groups <- 5
g.names <- c("Native","Asia","AUS-NZ","Europe","USA-CA")
g.codenames <- c("Native","Asia","AUS-NZ","Europe","USA-CA")
g.colors <- c("green", "blue", "orange","pink3", "red" )

# Import occurrences into data frames
setwd(main.directory)
occ.points.1 <- read.table("occurrences_native.txt", sep = "\t", header = TRUE)
occ.points.2 <- read.table("occurrences_asia.txt", sep = "\t", header = TRUE)
occ.points.3 <- read.table("occurrences_australia.txt", sep = "\t", header = TRUE)
occ.points.4 <- read.table("occurrences_europe.txt", sep = "\t", header = TRUE)
occ.points.5 <- read.table("occurrences_usa.txt", sep = "\t", header = TRUE)

# Join the occurrences into a list of data frames
occ.points <- list(occ.points.1,occ.points.2,occ.points.3,occ.points.4,occ.points.5)

# Import climatic variables information
variables.file.type <- '.bil'
variables.path <- "./variables"


# Offset that the minimum convex polygon (MCP) will apply to create the polygon
buffer.size <- 0.3

#Environmental space: An environmental space is generated based on the PCA values calculated for the background and the occurrence records. We defined the resolution of this two-dimensional space grid below
R <- 500

#Niche overlap: For the niche overlap, we calculate the D metric and its significance, using a similarity test. We define the number of interactions for the similarity test below (see the methods section in the manuscript for details)
rep <- 100

# ---------------- Preparation and visualization of the occurrence points  ----------------


# Thinning the occurrences. To reduce sampling bias, the occurrences are trimmed to ensure that a minimum distance of 5 km is lead between one to each other, so modify the parameter "thin.par" according to your study
occ.points.thin <- list()

for (i in 1:n.groups){
  occ.points.thin[i] <- thin(occ.points[[i]], verbose = FALSE, 
                          lat.col = "latitude",
                          long.col = "longitude",
                          spec.col = "species",
                          thin.par = 5,
                          reps = 1, 
                          write.files = FALSE,
                          write.log.file = FALSE,
                          locs.thinned.list.return = TRUE,
                          out.dir = getwd(),
                          out.base = "ThinProcess")
}


#To check the distribution of the occurrence records we map them in a world context
data(wrld_simpl)
plot(wrld_simpl)

for (i in 1:n.groups){
  points(occ.points.thin[[i]]$Longitude, occ.points.thin[[i]]$Latitude , col = g.colors[i], pch = 20, cex = 0.7)
}

#Next, we checked the number of occurrence records per region
number.occ <- list()

for (i in 1:n.groups) {
  number.occ[[i]] <- c(g.names[i],nrow(occ.points.thin[[i]]))
}

print(number.occ)  


# ---------------- Importation and Visualization of variables  ----------------

# Generate the paths vector to read
variables.file.pattern <- paste('*',variables.file.type, sep="")
list.variables<- list.files(path = variables.path, pattern = variables.file.pattern, full.names=TRUE)

# Import each variable file and store them in variables 
variables <- stack(list.variables)

# Plot of the variables
plot(variables)

# ---------------- Creation of polygons and their associated data (occurrences and background per region)  ----------------

# Define of the Minimum Convex Polygon (MCP) for each region based on their occurrences. This part of the code will take several minutes

# Function wich creates the MCP
mcp <- function (xy) {
  xy <- as.data.frame(coordinates(xy))
  coords.t <- chull(xy[, 1], xy[, 2])
  xy.bord <- xy[coords.t, ]
  xy.bord <- rbind(xy.bord[nrow(xy.bord), ], xy.bord)
  return(SpatialPolygons(list(Polygons(list(Polygon(as.matrix(xy.bord))), 1))))
}

# Union of the world map
lps <- getSpPPolygonsLabptSlots(wrld_simpl)
IDFourBins <- cut(lps[,1], range(lps[,1]), include.lowest=TRUE)
world <- unionSpatialPolygons(wrld_simpl, IDFourBins)

# Empty objects
xy.mcp <- list()   # List of Polygon data for each region
back.env <- list() # List of Background data for each region
spec.env <- list() # List of Occurrences data for each region

# Plot map
plot(wrld_simpl)

# Loop to create the MCP for each region

for(i in 1:n.groups) {

  # Background polygon creation
  mcp.occ <- mcp(occ.points.thin[[i]])                # Creation of the MCP without buffer
  xy.mcp.i <- gBuffer(mcp.occ, width = buffer.size)   # Application of the buffer size
  proj4string(xy.mcp.i) <- proj4string(world) 
  xy.mcp[[i]] <- gIntersection(xy.mcp.i, world, byid=TRUE, drop_lower_td=TRUE)

  # Background environment
  back.env[[i]] <- na.exclude(do.call(rbind.data.frame, extract(variables, xy.mcp[[i]])))

  # Species environment
  spec.env[[i]] <- na.exclude(extract(variables, occ.points.thin[[i]]))

  # Plot of the occurrences and the resultant MCP
  points(occ.points.thin[[i]], col = g.colors[i], pch = 20, cex = 0.7)
  plot(xy.mcp[[i]], add = TRUE, border = g.colors[i], lwd = 2)

  # Saving the MCP as a .shp (Shape file)
  xy.mcp[[i]] <- as(xy.mcp[[i]], "SpatialPolygonsDataFrame")

  shape.directory <- paste(getwd(),'/shapes',sep="")
  if (i==1){ dir.create(shape.directory) }
  setwd(shape.directory)
  shape.name <- paste("shape_", g.names[i], sep="")
  writeOGR(xy.mcp[[i]], dsn = '.', layer = shape.name, driver = "ESRI Shapefile")
  setwd(main.directory)
}


# ---------------- Preparation of the data to do the PCA  ----------------

# Environmental values for the background 
all.back.env <- do.call(rbind.data.frame, back.env)

# Environmental values for the species occurrence points 
all.spec.env <- do.call(rbind.data.frame, spec.env)

# Environmental values all together
data.env <- rbind(all.spec.env, all.back.env) 

# Saving the information of how all values are organized together in the data.env variable

# Number of rows and first row of the Species
n.rows.spec.env <- list()
first.row.spec.env <- list()

for (i in 1:n.groups){
  first.row.spec.env[i] <- do.call(sum, n.rows.spec.env)+1

  n.rows.spec.env[i] <- nrow(spec.env[[i]])
}

# Number of rows and first row of the Background
n.rows.back.env <- list()
first.row.back.env <- list()

for (i in 1:n.groups){
  first.row.back.env[i] <- do.call(sum, n.rows.spec.env)+1+do.call(sum,n.rows.back.env)
  n.rows.back.env[i] <- nrow(back.env[[i]])
} 

# Rows of the data.env that correspond to the Species
rows.spec.env <- list()

for (i in 1:n.groups){
    first.row <- first.row.spec.env[[i]]
    last.row <- first.row.spec.env[[i]] + n.rows.spec.env[[i]]-1
    rows.spec.env[[i]] <- c(first.row:last.row)
} 

# Rows of the data.env that correspond to the Background
rows.back.env <- list()

for (i in 1:n.groups){
  first.row <- first.row.back.env[[i]]
  last.row <- first.row.back.env[[i]] + n.rows.back.env[[i]]-1
  rows.back.env[[i]] <- c(first.row:last.row)
} 


# ---------------- PCA  ----------------

# Weight matrix ## 这个部分需要注意！
w <- c(rep(0, nrow(all.spec.env)), rep(1, nrow(all.back.env)))

# PCA of all environment
pca.cal <- dudi.pca(data.env, row.w = w, center = TRUE, 
                    scale = TRUE, scannf = FALSE, nf = 2)

# Empty lists to save the results
scores.spec <- list()
scores.back <- list()

# Assigning the results 
for(i in 1:n.groups) {
  scores.spec[[i]] <- pca.cal$li[rows.spec.env[[i]], ]
  scores.back[[i]] <- pca.cal$li[rows.back.env[[i]], ]  
}

# All the Background results together
total.scores.back <- do.call(rbind.data.frame, scores.back)


# ---------------- Species density in the grid  ----------------

#Next, we modeled the species density in the environmental grid, considering the observed occurrence density and the availability of the conditions in the background
z <- list()

for(i in 1:n.groups) {
  z[[i]] <- ecospat.grid.clim.dyn(total.scores.back,
                                  scores.back[[i]],
                                  scores.spec[[i]],
                                  R = R)
}

# ---------------- Niche overlap  ---------------- ####

#Niche overlap: For the niche overlap, we calculate the D metric and its significance, using a similarity test. The number of iterations has been defined at the beginning of the code

#Once the number of interactions is defined, we can generate the simulated overlap D values. Additionally, we calculate the niche dynamic indices: niche unfilling, expansion, and stability (see methods in the manuscript)

# Empty matrices
D <- matrix(nrow = n.groups, ncol = n.groups)
rownames(D) <- colnames(D) <- g.codenames
unfilling <- stability <- expansion <- sim.d <- sim.s <- eq.d <- eq.s <- D

# Filling the matrices
for(i in 2:n.groups) {

  for(j in 1:(i - 1)) {

    x1 <- z[[i]]
    x2 <- z[[j]]

    # Niche overlap
    D[i, j] <- ecospat.niche.overlap (x1, x2, cor = TRUE)$D

    # Niche equivalency. Note: The quantile of the environmental density used to remove marginal climates. If             intersection=NA, the analysis is performed on the whole environmental extent (native and invaded). If intersection     =0, the analysis is performed at the intersection between native and invaded range. If intersection=0.05, the         analysis is performed at the intersection of the 5th quantile of both native and invade environmental densitie
    eq.s[i, j] <- ecospat.niche.equivalency.test (x1, x2, rep, alternative = "greater")$p.I
    eq.s[j, i] <- ecospat.niche.equivalency.test (x2, x1, rep, alternative = "greater")$p.I

    eq.d[i, j] <- ecospat.niche.equivalency.test (x1, x2, rep, alternative = "lower")$p.I
    eq.d[j, i] <- ecospat.niche.equivalency.test (x2, x1, rep, alternative = "lower")$p.I


    # Niche similarity to detect more similar niches than expected by change (i.e. niche overlap is more equivalent       /similar than random)
    sim.s[i, j] <- ecospat.niche.similarity.test (x1, x2, rep, alternative = "greater")$p.D
    sim.s[j, i] <- ecospat.niche.similarity.test (x2, x1, rep, alternative = "greater")$p.D

    # Niche similarity to detect more different niches than expected by chance (i.e. the niche overlap is less            equivalent/similar than random)
    sim.d[i, j] <- ecospat.niche.similarity.test (x1, x2, rep, alternative = "lower")$p.D
    sim.d[j, i] <- ecospat.niche.similarity.test (x2, x1, rep, alternative = "lower")$p.D


    # Niche Expansion, Stability, and Unfilling
    index1 <- ecospat.niche.dyn.index (x1, x2, 
                                       intersection = NA)$dynamic.index.w
    index2 <- ecospat.niche.dyn.index (x2, x1,
                                       intersection = NA)$dynamic.index.w
    expansion[i, j] <- index1[1]
    stability[i, j] <- index1[2]
    unfilling[i, j] <- index1[3]
    expansion[j, i] <- index2[1]
    stability[j, i] <- index2[2]
    unfilling[j, i] <- index2[3]
  }
}

#Numeric results: Below, we present the complete results for each metric, considering all pair-wise comparisons between all range areas

#D value
kable(D, digits = 3, format = "markdown")

#Niche equivalency null model (p-values) to detect equivalent niches
kable(eq.s, digits = 3, format = "markdown")

#Niche equivalency null model (p-values) to detect non equivalent niches
kable(eq.d, digits = 3, format = "markdown")

#Niche similarity null model (p-values) to detect similar niches
kable(sim.s, digits = 3, format = "markdown")

#Niche similarity null model (p-values) to detect different niches
kable(sim.d, digits = 3, format = "markdown")

#Niche Unfilling:
kable(unfilling, digits = 3,  format = "markdown")

#Niche Expansion:
kable(expansion, digits = 3,  format = "markdown")

#Niche Stability:
kable(stability, digits = 3,  format = "markdown")



# ---------------- Plots  ----------------####

# Function to plot the data
plot.niche.all <- function(z, n.groups, g.names,

                           densidade,
                           n.degradation,
                           slope.degradation,

                           contornar,
                           prob.dens,
                           line.width.dens,

                           back,
                           prob.back.1,
                           prob.back.2,
                           line.width.back,

                           title,
                           g.colors, 
                           add.plots, i) {  


  color.vector <- function(g.colors,n.degradation,prob.dens){

    alpha.seq <- numeric(length=n.degradation+1)
    alpha.x <- seq(0,1,(1/n.degradation))

    for (j in 1:(n.degradation+1)){
      alpha.seq[j] <- (slope.degradation^alpha.x[j]-1)/(slope.degradation-1)
    }
    col.vector <- numeric(n.degradation+1)
    current.color <- col2rgb(g.colors)/255
    for (j in 1:(n.degradation+1)){
      col.vector[j] <- rgb(current.color[1, ], current.color[2, ], current.color[3, ], alpha = alpha.seq[j])
    }
    return(col.vector)
  }

  a <- colorRampPalette( c(color.vector(g.colors[i], n.degradation)), alpha = TRUE )  


  xlim <- c(min(sapply(z, function(x){min(x$x)})),
            max(sapply(z, function(x){max(x$x)})))

  ylim <- c(min(sapply(z, function(x){min(x$y)})),
            max(sapply(z, function(x){max(x$y)})))

  z[[i]]$z.uncor@data@values <- matrix(data = z[[i]]$z.uncor@data@values, nrow = R, ncol = R, byrow = TRUE)
  z[[i]]$Z@data@values <- matrix(data = z[[i]]$Z@data@values, nrow = R, ncol = R, byrow = TRUE)

  image(z[[i]]$x, z[[i]]$y, z[[i]]$z.uncor@data@values, col = "transparent", 
        ylim = ylim, xlim = xlim,
        zlim = c(0.000001, max(z[[1]]$Z@data@values, na.rm = T)), 
        xlab = "PC1", ylab = "PC2",cex.lab = 1.5,
        cex.axis = 1.4,
        add = add.plots)
  abline(h = 0, v = 0, lty = 2)
  box()
  title(title)

  if (back) {
    contour(z[[i]]$x, z[[i]]$y, z[[i]]$Z@data@values, add = back,
            levels = quantile(z[[i]]$Z@data@values[z[[i]]$Z@data@values > 0], probs = 1-prob.back.1),
            drawlabels = FALSE,lty =1,
            col = g.colors[i], lwd = line.width.back)
    contour(z[[i]]$x, z[[i]]$y, z[[i]]$Z@data@values, add = back,
            levels = quantile(z[[i]]$Z@data@values[z[[i]]$Z@data@values > 0], probs = 1-prob.back.2),
            drawlabels = FALSE,lty =2,
            col = g.colors[i], lwd = line.width.back)
  }

  if (densidade) {

    densidade.value <- quantile(z[[i]]$z.uncor@data@values[z[[i]]$z.uncor@data@values > 0],probs = 1-prob.dens)

    densidade.quantile <- z[[i]]$z.uncor@data@values

    densidade.quantile[densidade.quantile < densidade.value] <- 0

    densidade.quantile <- matrix(data = densidade.quantile, nrow = R, ncol = R, byrow = FALSE)

    image(z[[i]]$x, z[[i]]$y, densidade.quantile, col = a(n.degradation+1), add = densidade)

  }


  if(contornar){
    contour(z[[i]]$x, z[[i]]$y, z[[i]]$z.uncor@data@values, 
            add = contornar,
            levels = quantile(z[[i]]$z.uncor@data@values[z[[i]]$z.uncor@data@values > 0],probs = 1-prob.dens),
            drawlabels = FALSE,
            col = g.colors[i],
            lwd = line.width.dens, lty = 1)
  }
}


# Individual niche plots

for(i in 1:n.groups) {
  plot.niche.all(z, n.groups, g.names,

                 densidade = TRUE,          # TRUE/FALSE to activate/desactivate the density plot
                 n.degradation = 1000,      # Color resolution
                 slope.degradation = 100,   # Color degradation constant. It has to be >1

                 contornar = FALSE,         # TRUE/FALSE to activate/desactivate the contour of %density plot
                 prob.dens = 0.5,           # % of the density that represents the contour (0->0% and 1->100%)
                 line.width.dens = 2,       # Line width of the contour of %density

                 back = TRUE,               # TRUE/FALSE to activate/desactivate the background plot
                 prob.back.1 = 1,           # % of the background that represents the contour 1 (0->0% and 1->100%)
                 prob.back.2 = 0.5,         # % of the background that represents the contour 2 (0->0% and 1->100%)
                 line.width.back = 2,       # Line width of the contour of background

                 title = g.names[i],        # Title of each region
                 g.colors,                  # Colors of the regions
                 add.plots = FALSE, i)      # TRUE/FALSE to activate/desactivate the adition of plots in the same image
}

# Plot global Contour %Density + Occurrences Densities

for(i in 1:n.groups) {
  if (i==1){ add.plots <- FALSE }
  else { add.plots <- TRUE}

  plot.niche.all(z, n.groups, g.names,

                 densidade = TRUE,           # TRUE/FALSE to activate/desactivate the density plot
                 n.degradation = 1000,       # Color resolution
                 slope.degradation = 100,    # Color degradation constant. It has to be >1

                 contornar = TRUE,           # TRUE/FALSE to activate/desactivate the contour of %density plot
                 prob.dens = 0.2,            # % of the density that represents the contour (0->0% and 1->100%)
                 line.width.dens = 2,        # Line width of the contour of %density

                 back = FALSE,               # TRUE/FALSE to activate/desactivate the background plot
                 prob.back.1 = 1,           # % of the background that represents the contour 1 (0->0% and 1->100%)
                 prob.back.2 = 0.5,         # % of the background that represents the contour 2 (0->0% and 1->100%)
                 line.width.back = 2,        # Line width of the contour of background

                 title = "Densities",        # Title of the plot
                 g.colors,                   # Colors of the regions
                 add.plots = add.plots, i)   # TRUE/FALSE to activate/desactivate the adition of plots in the same image
}

# Plot global Contour %Density + Occurrences Densities + Background climate

for(i in 1:n.groups) {
  if (i==1){ add.plots <- FALSE }
  else { add.plots <- TRUE}

  plot.niche.all(z, n.groups, g.names,

                 densidade = TRUE,          # TRUE/FALSE to activate/desactivate the density plot
                 n.degradation = 1000,      # Color resolution
                 slope.degradation = 100,   # Color degradation constant. It has to be >1

                 contornar = TRUE,          # TRUE/FALSE to activate/desactivate the contour of %density plot
                 prob.dens = 0.98,           # % of the density that represents the contour (0->0% and 1->100%)
                 line.width.dens = 1,       # Line width of the contour of %density

                 back = TRUE,               # TRUE/FALSE to activate/desactivate the background plot
                 prob.back.1 = 0.98,           # % of the background that represents the contour 1 (0->0% and 1->100%)
                 prob.back.2 = FALSE,       # % of the background that represents the contour 2 (0->0% and 1->100%)
                 line.width.back = 3,       # Line width of the contour of background

                 title = "Densities and Background",        # Title of the plot
                 g.colors,                   # Colors of the regions
                 add.plots = add.plots, i)   # TRUE/FALSE to activate/desactivate the adition of plots in the same image 
}


# Loadings plot ####

#Below the loadings plot (contribution of the variables for each axis). Check the variable codes at http://www.worldclim.org/bioclim

loadings <- cbind(cor(data.env, pca.cal$tab[,1]), cor(data.env, pca.cal$tab[,2]))
colnames(loadings) <- c("axis1", "axis2")

# Plot of the contribution of each variable to the PC1
barplot(loadings[,1], las=2, main="PC1")

# Plot of the contribution of each variable to the PC2
barplot(loadings[,2], las=2, main="PC2")


# Arrows plot

# Plot of the arrows representing the contribution of each variable, directly on the environmental space

# Get data from PCA
contrib <- pca.cal$co
eigen <- pca.cal$eig

# Preparation of the names of each variable
nomes <- gsub(paste(variables.path,'/', sep=""),'',list.variables)
nomes <- gsub(variables.file.type,'',nomes)

# Plot of the circle
s.corcircle(contrib[, 1:2] / max(abs(contrib[, 1:2])), grid = F,  label = nomes, clabel = 1.2)

# Addition of the axis information
text(0, -1.1, paste("PC1 (", round(eigen[1]/sum(eigen)*100,2),"%)", sep = ""))
text(1.1, 0, paste("PC2 (", round(eigen[2]/sum(eigen)*100,2),"%)", sep = ""), srt = 90)